%% GranularResidualCES_bycountry.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Date: May 4, 2019
%
% Program: this is a first attempt at operationalizing the granular
% residual calculations following from master_notesv3.pdf (page 34) using
% the methodology of Gabaix and Koijen's "Granular IV" paper. To begin, we
% will do this based on PWT data for real GDP growth as well as TFP
% measures reported in this dataset. 
%
% The PWT dataset will be for 1970-2014, and will be loaded as levels, and
% then manipulated.
%
% Once the granular residuals are constructed based on the different
% partner country growth rates, we will compare to France's actual real GDP
% and TFP growth to see if any patterns exist. I.e., the model generated
% elasticities that can map the country shocks into an aggregate foreign
% component that can predct overall French GDP or TFP growth.
%
% We can also compare to the granular residuals we've calculated from
% previous work using French firm-level data, but will need to restrict the
% sample size over years further for this.
%
% UPDATED 2/21/20 (JDG): adjusted calls of directories given still in
% Robustness folder. Will need to be adjustd later. Also switched to double
% deflation method for calculating elasticities.
%
% Also created a version of Table 4 of paper to output.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
close all

%%
% Loan PWT real gdp, exchange rate, and TFP data from 1970-2014

% rgdpe   Expenditure-side real GDP at chained PPPs (in mil. 2011US$)
% rgdpo   Output-side real GDP at chained PPPs (in mil. 2011US$)
% rgdpna  Real GDP at constant 2011 national prices (in mil. 2011US$)
% ctfp    TFP level at current PPPs (USA=1)
% cwtfp   Welfare-relevant TFP levels at current PPPs (USA=1)
% xr      Exchange rate, national currency/USD (market+estimated)
% 
 load ../../../Data/other/pwt9_rgdpftp.mat
% 
 ISO = string(pwt9rgdptfp.ISO);
 year_pwt = pwt9rgdptfp.year;
% rgdpe = pwt9rgdptfp.rgdpe;
% rgdpo = pwt9rgdptfp.rgdpo;
% rgdpna = pwt9rgdptfp.rgdpna;
 ctfp = pwt9rgdptfp.ctfp;
 cwtfp = pwt9rgdptfp.cwtfp;
% xr = pwt9rgdptfp.xr;


load ../../../Data/RGDPlcu_WDI.mat

ISO_WDI = string(RGDPlcuWDIclean.ISO);
year_WDI = RGDPlcuWDIclean.year;
rgdpe = RGDPlcuWDIclean.Value;


%%
% Caculate growth rates of rgdpe (in local currency_ and ctfp
% Will fill in NaN for missing values and keep a balanced panel for now

% Transform RGDP to local currency to get rid of exchange rate swings for

%rgdpe_lc = rgdpe.*xr;
rgdpe_lc = rgdpe;
% A numerical index to loop over countries

id = grp2idx(ISO);
id_WDI = grp2idx(ISO_WDI);

france = 15;

% Initialize vector of growth rates for real GDP and TFP

gr_rgdpe_lc = NaN(size(id_WDI,1),1);
gr_ctfp = NaN(size(id,1),1);

for i = 1:max(id)
   
   x = ctfp(id==i);
   grx = NaN(size(x,1),1);
   grx(2:end) = x(2:end)./x(1:end-1)-1;
   gr_ctfp(id==i) = grx;   
   
   clear x grx;
end

for i = 1:max(id_WDI)
   x = rgdpe_lc(id_WDI==i);
   grx = NaN(size(x,1),1);
   grx(2:end) = x(2:end)./x(1:end-1)-1;
   gr_rgdpe_lc(id_WDI==i) = grx;
   
   clear x grx;
   
end

%%
% Calculate Granular residual based on partner countries:
%
%  a) TFP growth
%  b) Real GDP growth
%
% Note: the numnber of observations may vary per country, and we
% will compare to overall GDP and TFP growth of France later.

CTY = {'AUS', 'AUT', 'BEL', 'BGR', 'BRA', 'CAN', 'CHN' ,'CYP', ...
'CZE', 'DEU', 'DNK', 'ESP' , 'EST', 'FIN', 'FRA', 'GBR', ...
'GRC', 'HUN', 'IDN', 'IND', 'IRL', 'ITA', 'JPN', 'KOR', ...
'LTU' ,'LVA', 'MEX','MLT', 'NLD', 'POL', 'PRT', 'ROU','RUS', ...
'SVK', 'SVN', 'SWE', 'TUR', 'TWN', 'USA','ROW'};
% Initialize vectors we will collect results in

T = size(ctfp,1)/max(id);
T_WDI = size(rgdpe_lc,1)/max(id_WDI);
Gamma_nRGDP = zeros(T_WDI,size(CTY,2));
Gamma_nTFP = zeros(T,size(CTY,2));
Overall_nRGDP = zeros(T_WDI,size(CTY,2));
Overall_nTFP = zeros(T,size(CTY,2));

% Run loop a calculate Granular Residula w.r.t. to each country's shock. 
% NOTE: only do w.r.t. TFP growth now given we need to adjust formula for
% real GDP growth calculation -- imperfect what we do for GDP

for i = 1:size(CTY,2)
   cty = string(CTY(i));
   if cty~='FRA' & cty~='ROW'
   % Get required counterfactual data
   
   fname = strcat('MAT/rho3eta1.000001lambda1.000001_psi3/Step_5_output_alpha4_rho3_',cty,'prod_p10_approx.mat');   
   load(fname,'VA_fj_0','VA_hat_fj','RGDP_cpi_hat_n','RGDP_def_hat_n','GDP_deflator_n','DoubleDeflation_deflator_n','RGDP_doubledef_hat_n');   
   
   % Calculate shares for weighted and unweighted aggregation
   
   w_sh = VA_fj_0/sum(VA_fj_0);
   uw_sh = 1/size(VA_fj_0,1);
   
   % Calculate elastitcity at firm level given 10% productivity shock to
   % country in counterfactual
   
   e_f = (VA_hat_fj-1)/.1;          % F*1
   e_f_real = (VA_hat_fj/DoubleDeflation_deflator_n(france)-1)/.1;          % F*1
   e_f_GDP = (VA_hat_fj/DoubleDeflation_deflator_n(france)-1)./(RGDP_doubledef_hat_n(i,1)-1); 
%    e_f_real = (VA_hat_fj/GDP_deflator_n(france)-1)/.1;          % F*1
%    e_f_GDP = (VA_hat_fj/GDP_deflator_n(france)-1)./(RGDP_def_hat_n(i,1)-1);
   cty
   mean(e_f)
   mean(e_f_real)
   mean(e_f_GDP)
   % Calculate Granular Residual component per country per time period. 
   % Given we are summing over firms and countries, we will sum over firms
   % per country in loop and then we will sum across countries after loop
   % to get total Granular Residual

   % (a) TFP Growth Granular Residual
   
   growth =  gr_ctfp(ISO==cty);      % T*1
   
   for t = 1:T
      granres_n = sum(w_sh.*e_f_real*growth(t)) - mean(e_f_real*growth(t));
      Gamma_nTFP(t,i) = granres_n;
      Overall_nTFP(t,i) = sum(w_sh.*e_f_real*growth(t));
      Mean_nTFP(t,i) = mean(e_f_real*growth(t));
      Mean_nTFP_(t,i) = Overall_nTFP(t,i)-Gamma_nTFP(t,i);
      clear granres_n;
   end
   
   clear growth
   
   % (b) Real GDP Growth Granular Residual
   if cty~='TWN'
   growth =  gr_rgdpe_lc(ISO_WDI==cty);      % T*1
   
   for t = 1:T_WDI
      granres_n = sum(w_sh.*e_f_GDP*growth(t)) - mean(e_f_GDP*growth(t));
      Gamma_nRGDP(t,i) = granres_n;
      Overall_nRGDP(t,i) = sum(w_sh.*e_f_GDP*growth(t));
      Mean_nRGDP(t,i) = mean(e_f_GDP*growth(t));
      Mean_nRGDP_(t,i) = Overall_nRGDP(t,i)-Gamma_nRGDP(t,i);
      clear granres_n;
   end
   end
   clear cty fname w_sh uw_sh e_f growth
   end
end

year_WDI_=unique(year_WDI);
year_pwt_=unique(year_pwt);

Gamma_TFP = nansum(Gamma_nTFP,2);
Gamma_RGDP = nansum(Gamma_nRGDP,2);
Overall_TFP = nansum(Overall_nTFP,2);
Overall_RGDP = nansum(Overall_nRGDP,2);
Mean_TFP = nansum(Mean_nTFP,2);
Mean_RGDP = nansum(Mean_nRGDP,2);
Mean_TFP_ = nansum(Mean_nTFP_,2);
Mean_RGDP_ = nansum(Mean_nRGDP_,2);

std(Overall_TFP(2:end))
std(Mean_TFP(2:end))
std(Mean_TFP_(2:end))
std(Gamma_TFP(2:end))

grTFP_FRA = gr_ctfp(ISO=="FRA");
grRGDP_FRA = gr_rgdpe_lc(ISO_WDI=="FRA");

% Summary statistics and figures
ratio_Gamma_TFP=Gamma_TFP(12:end)./grRGDP_FRA(2:55);
ratio_Gamma_RGDP=Gamma_RGDP./grRGDP_FRA;
nanmedian(ratio_Gamma_TFP)
nanmedian(ratio_Gamma_RGDP)
nanmean(ratio_Gamma_TFP)
nanmean(ratio_Gamma_RGDP)


b = regstats(Gamma_TFP(12:end),grRGDP_FRA(2:55),'linear','yhat');
figure(1)
plot(grRGDP_FRA(2:55),Gamma_TFP(12:end),'bo',grRGDP_FRA(2:55),b.yhat,'r-','LineWidth',2)
xlabel('GDP Growth')
ylabel('TFP Granular Residual')
saveas(gcf,'Figures/rho3eta1.000001lambda1.000001_psi3/Gran/TFPGran_vs_GDPGrowth.pdf')
clear b
std(Gamma_TFP(12:end))
std(grRGDP_FRA(2:55))
corr(Gamma_TFP(12:end),grRGDP_FRA(2:55))
std(Gamma_TFP(42:58))
std(grRGDP_FRA(32:48))
corr(Gamma_TFP(42:58),grRGDP_FRA(32:48))

b = regstats(Gamma_RGDP(2:end-1),grRGDP_FRA(2:end-1),'linear','yhat');
figure(2)
plot(grRGDP_FRA(2:end-1),Gamma_RGDP(2:end-1),'bo',grRGDP_FRA(2:end-1),b.yhat,'r-','LineWidth',2)
xlabel('GDP Growth')
ylabel('GDP Granular Residual')
saveas(gcf,'Figures/rho3eta1.000001lambda1.000001_psi3/Gran/GDPGran_vs_GDPGrowth.pdf')
clear b
std(Gamma_RGDP(2:end-1))
std(grRGDP_FRA(2:end-1))
corr(Gamma_RGDP(2:end-1),grRGDP_FRA(2:end-1))
std(Gamma_RGDP(15:end-1))
std(grRGDP_FRA(15:end-1))
corr(Gamma_RGDP(15:end-1),grRGDP_FRA(15:end-1))

yyaxis left
plot(year_WDI_(2:end-1),Gamma_RGDP(2:end-1),'r-');
ylabel('GDP Granular Residual')
yyaxis right
plot(year_WDI_(2:end-1),grRGDP_FRA(2:end-1),'b-');
ylabel('GDP Growth')
xlabel('')
saveas(gcf,'Figures/rho3eta1.000001lambda1.000001_psi3/Gran/GDPGran_GDPGrowth_line.pdf')


yyaxis left
plot(year_WDI_(2:55),Gamma_TFP(12:end),'r-');
ylabel('TFP Granular Residual')
yyaxis right
plot(year_WDI_(2:55),grRGDP_FRA(2:55),'b-');
ylabel('GDP Growth')
xlabel('')
saveas(gcf,'Figures/rho3eta1.000001lambda1.000001_psi3/Gran/TFPGran_GDPGrowth_line.pdf')


% Compute granular residual for France based on results of the Econometrica
% paper

load('MAT/reg_output_all_ca');
wgr_idio = w_fp.*gr_idio;
[agr_idio sd_agidio sgr_idio] = ag_growth(id_fp,year,wgr_idio,gr_idio);
agr_idio=agr_idio';
year_ficus=unique(year);

b = regstats(Gamma_TFP(42:58),agr_idio,'linear','yhat');
labels = cellstr(num2str(year_ficus));
figure(1)
plot(agr_idio,Gamma_TFP(42:58),'bo',agr_idio,b.yhat,'r-','LineWidth',2)
%text(agr_idio,Gamma_TFP(42:58),labels,'VerticalAlignment','bottom','HorizontalAlignment','right')
xlabel('French Granular Residual')
ylabel('TFP Granular Residual')
saveas(gcf,'Figures/rho3eta1.000001lambda1.000001_psi3/Gran/TFPGran_vs_FrenchGran.pdf')
clear b
std(grRGDP_FRA(32:48))
std(agr_idio)
std(Gamma_TFP(42:58))
corr(Gamma_TFP(42:58),grRGDP_FRA(32:48))
corr(Gamma_TFP(42:58),agr_idio)

b = regstats(Gamma_RGDP(32:48),agr_idio,'linear','yhat');
figure(2)
plot(agr_idio,Gamma_RGDP(32:48),'bo',agr_idio,b.yhat,'r-','LineWidth',2)
%text(agr_idio,Gamma_RGDP(32:48),labels,'VerticalAlignment','bottom','HorizontalAlignment','right')
xlabel('French Granular Residual')
ylabel('GDP Granular Residual')
saveas(gcf,'Figures/rho3eta1.000001lambda1.000001_psi3/Gran/GDPGran_vs_FrenchGran.pdf')
clear b
std(grRGDP_FRA(32:48))
std(agr_idio)
std(Gamma_RGDP(32:48))
corr(Gamma_RGDP(32:48),grRGDP_FRA(32:48))
corr(Gamma_RGDP(32:48),agr_idio)

%% Generate table of output using write table

% Data moments

sdRGDPgr_7514 = std(grRGDP_FRA(year_WDI_>=1975&year_WDI_<=2014))*100;
sdRGDPgr_9107 = std(grRGDP_FRA(year_WDI_>=1991&year_WDI_<=2007))*100;
stGRANres_9107 = sd_agidio*100;

% Foreign TFP moments

sddlYF_TFP_7514 = std(Overall_TFP(year_pwt_>=1975&year_pwt_<=2014))*100;
sddGF_TFP_7514 = std(Gamma_TFP(year_pwt_>=1975&year_pwt_<=2014))*100;
sdEF_TFP_7514 = std(Mean_TFP(year_pwt_>=1975&year_pwt_<=2014))*100;

sddlYF_TFP_9107 = std(Overall_TFP(year_pwt_>=1991&year_pwt_<=2007))*100;
sddGF_TFP_9107 = std(Gamma_TFP(year_pwt_>=1991&year_pwt_<=2007))*100;
sdEF_TFP_9107 = std(Mean_TFP(year_pwt_>=1991&year_pwt_<=2007))*100;

% Foreign GDP moments

sddlYF_RGDP_7514 = std(Overall_RGDP(year_WDI_>=1975&year_WDI_<=2014))*100;
sddGF_RGDP_7514 = std(Gamma_RGDP(year_WDI_>=1975&year_WDI_<=2014))*100;
sdEF_RGDP_7514 = std(Mean_RGDP(year_WDI_>=1975&year_WDI_<=2014))*100;

sddlYF_RGDP_9107 = std(Overall_RGDP(year_WDI_>=1991&year_WDI_<=2007))*100;
sddGF_RGDP_9107 = std(Gamma_RGDP(year_WDI_>=1991&year_WDI_<=2007))*100;
sdEF_RGDP_9107 = std(Mean_RGDP(year_WDI_>=1991&year_WDI_<=2007))*100;

% Collect terms for table

table_data = [sdRGDPgr_7514 NaN;sdRGDPgr_9107 stGRANres_9107];
table_tfp = [sddlYF_TFP_7514 sddGF_TFP_7514 sdEF_TFP_7514;sddlYF_TFP_9107 sddGF_TFP_9107 sdEF_TFP_9107];
table_rgdp = [sddlYF_RGDP_7514 sddGF_RGDP_7514 sdEF_RGDP_7514; sddlYF_RGDP_9107 sddGF_RGDP_9107 sdEF_RGDP_9107];

table_granres = table(table_data,table_tfp,table_rgdp);
table_granres.Properties.RowNames = {'1975--2014' '1991--2007'};

writetable(table_granres,'Tables/table_granres.xlsx','WriteVariableNames',true,'WriteRowNames',true);
% Foreign GDP moments

% Use granular residual for France based on our panel of firm-level data

%
% DATA = dlmread('CSV/Granular_2.csv',',');
% year_panel=DATA(:,1);
% Gamma_vabcf=DATA(:,2);
% agr_idio_ficus=DATA(:,3);
% 
% b = regstats(Gamma_TFP(47:58),agr_idio_ficus,'linear','yhat');
% labels = cellstr(num2str(year_panel));
% figure(1)
% plot(agr_idio_ficus,Gamma_TFP(47:58),'bo',agr_idio_ficus,b.yhat,'r-','LineWidth',2)
% text(agr_idio_ficus,Gamma_TFP(47:58),labels,'VerticalAlignment','bottom','HorizontalAlignment','right')
% xlabel('French Granular Residual')
% ylabel('TFP Granular Residual')
% clear b
% corr(Gamma_TFP(47:58),agr_idio_ficus)
% 
% b = regstats(Gamma_TFP(47:58),Gamma_vabcf,'linear','yhat');
% labels = cellstr(num2str(year_panel));
% figure(1)
% plot(Gamma_vabcf,Gamma_TFP(47:58),'bo',Gamma_vabcf,b.yhat,'r-','LineWidth',2)
% text(Gamma_vabcf,Gamma_TFP(47:58),labels,'VerticalAlignment','bottom','HorizontalAlignment','right')
% xlabel('French Granular Residual')
% ylabel('TFP Granular Residual')
% clear b
% corr(Gamma_TFP(47:58),Gamma_vabcf)
% 
% b = regstats(Gamma_RGDP(37:48),agr_idio_ficus,'linear','yhat');
% labels = cellstr(num2str(year_panel));
% figure(1)
% plot(agr_idio_ficus,Gamma_RGDP(37:48),'bo',agr_idio_ficus,b.yhat,'r-','LineWidth',2)
% text(agr_idio_ficus,Gamma_RGDP(37:48),labels,'VerticalAlignment','bottom','HorizontalAlignment','right')
% xlabel('French Granular Residual')
% ylabel('GDP Granular Residual')
% clear b
% corr(Gamma_RGDP(37:48),agr_idio_ficus)
% 
% b = regstats(Gamma_RGDP(37:48),Gamma_vabcf,'linear','yhat');
% labels = cellstr(num2str(year_panel));
% figure(1)
% plot(Gamma_vabcf,Gamma_RGDP(37:48),'bo',Gamma_vabcf,b.yhat,'r-','LineWidth',2)
% text(Gamma_vabcf,Gamma_RGDP(37:48),labels,'VerticalAlignment','bottom','HorizontalAlignment','right')
% xlabel('French Granular Residual')
% ylabel('GDP Granular Residual')
% clear b
% corr(Gamma_RGDP(37:48),Gamma_vabcf)



% 
% plot(year_WDI_(2:end-1),Gamma_RGDP(2:end-1),'r-',year_WDI_(2:end-1),grRGDP_FRA(2:end-1),'g-')
% % Plot growth of French GDP and compare with INSEE data 
% % https://www.insee.fr/fr/statistiques/serie/010548503#Graphique
% load ../../Data/INSEE_rgdp.mat
% grRGDP_FRA_INSEE=INSEE_RGDPp(:,2)/100;
% year_=unique(year_pwt);
% grRGDP_FRA_comp=[grRGDP_FRA grRGDP_FRA_INSEE];
% % Create figure
% figure1 = figure;
% % Create axes
% axes1 = axes('Parent',figure1);
% hold(axes1,'on');
% 
% % Create multiple lines using matrix input to plot
% plot1 = plot(year_,grRGDP_FRA_comp,'LineWidth',2);
% set(plot1(1),'DisplayName','PWT','Color',[1 0 0]);
% set(plot1(2),'DisplayName','INSEE','Color',[0 0 1]);
% 
% % Create ylabel
% ylabel({'Real GDP growth'});
% 
% % Uncomment the following line to preserve the X-limits of the axes
% % xlim(axes1,[1950 2020]);
% % Uncomment the following line to preserve the Y-limits of the axes
% % ylim(axes1,[-0.05 0.1]);
% box(axes1,'on');
% % Create legend
% legend1 = legend(axes1,'show');
% 
% plot(year_(2:end),grRGDP_FRA(2:end,:),'r-','Linewidth',2);
% hold on;
% plot(year_(2:end),grRGDP_FRA_INSEE(2:end,:),'b-','Linewidth',2);
% hold off;
% 
